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The problem of the detection and mapping of a stochastic gravitational wave background (SGWB), 
either cosmological or astrophysical, bears a strong semblance to the analysis of the cosmic microwave 
background (CMB) anisotropy and polarization, which too is a stochastic field, statistically described 
in terms of its correlation properties. An astrophysical gravitational wave background (AGWB) will 
likely arise from an incoherent superposition of unmodelled and/or unresolved sources and cosmo- 
logical gravitational wave backgrounds (CGWB) are also predicted in certain scenarios. The basic 
statistic we use is the cross-correlation between the data from a pair of detectors. In order to 'point' 
the pair of detectors at different locations one must suitably delay the signal by the amount it takes 
for the gravitational waves (GW) to travel to both detectors corresponding to a source direction. 
Then the raw (observed) sky map of the SGWB is the signal convolved with a beam response func- 
tion that varies with location in the sky. We first present a thorough analytic understanding of the 
structure of the beam response function using an analytic approach employing the stationary phase 
approximation. The true sky map is obtained by numerically deconvolving the beam function in 
the integral (convolution) equation. We adopt the maximum likelihood framework to estimate the 
true sky map using the conjugate gradient method that has been successfully used in the broadly 
similar, well-studied CMB map making problem. We numerically implement and demonstrate the 
method on signal generated by simulated (unpolarized) SGWB for the GW radiometer consisting 
of the LIGO pair of detectors at Hanford and Livingston. We include 'realistic' additive Gaussian 
noise in each data stream based on the LIGO-I noise power spectral density. The extension of the 
method to multiple baselines and polarized GWB is outlined. In the near future the network of GW 
detectors, including the Advanced LIGO and Virgo detectors that will be sensitive to sources within 
a thousand times larger spatial volume, could provide promising data sets for GW radiometry. 

PACS numbers: 04.80.Nn, 04.30. Db, 95.55.Ym, 98.70.Vc, 98.80. Es, 07.05.Kf, 95.75.Pq 



I. INTRODUCTION 

The existence of gravitational waves (GW), has long 
been verified 'indirectly' through the observations of 
Hulse and Taylor [1]. However, direct observation of such 
waves with manmade gravitational wave detectors has 
been lacking. At present the laser interferometric detec- 
tors have achieved sensitivities close to that required for 
detecting such waves [2]. The space mission LISA [3] is 
also planned by the ESA and NASA to detect low fre- 
quency GW. The significance of the direct detection of 
GW lies, not only in the opening of an entirely new win- 
dow into observational astronomy by probing phenomena 
in the regime of strong gravity; it further promises to test 
our present theories of gravitation. 
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Different types of GW sources have been predicted and 
may be directly observed by Earth-based detectors in the 
near future (see [4, 5, 6, 7, 8, 9, 10] and references therein 
for recent reviews): (i) Transient sources - such as binary 
systems of neutron stars (NS) and/or black holes (BH) 
in their in-spiral phase, BH/BH and/or BH/NS merg- 
ers, and supernovae explosions, whose signals last for a 
time much shorter, typically between a few milli-seconds 
and a few minutes, than the planned observational time; 
(ii) Continuous wave (CW) sources - e.g. rapidly rotat- 
ing neutron stars, where a weak deterministic signal is 
continuously emitted, and (iii) Stochastic backgrounds 
of radiation, either of primordial or astrophysical origin. 

In this paper we will address the problem of a spa- 
tially resolved search of the gravitational wave stochastic 
background. This approach was advocated in the LIGO 
technical note [11] and the basic analysis was recently 
implemented on the fourth science run data from the 
LIGO interferometers to prepare an upper limit map [12]. 
Our main focus will be on a stochastic astrophysical GW 
background (AGWB), which might arise from a superpo- 
sition of a large number of independent and unresolved 
GW sources. The gravitational wave background can 
arise from a variety of sources: supernovae with asym- 
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metric core collapse, binary black hole (BBH) mergers, 
GWs from low-mass X-ray binaries (LMXBs) and hy- 
drodynamical instabilities in neutron stars (r-modes), 
or even GWs from astrophysical objects that we never 
knew existed. When a collection of any subset of these 
sources is unresolvable, it can appear as a stochastic GW 
background (SGWB) of a variable duration in our detec- 
tors of interest. While an astrophysical background will 
provide information about our immediate neighborhood, 
cosmological GW backgrounds (CGWB) could probe the 
physics of the early universe. There exist cosmological 
scenarios (e.g., cosmic strings and super-string models) 
which predict CGWB that should be detectable by Ad- 
vanced LIGO [13]. 

We propose and develop a data analysis method that 
measures and maps the power in the SGWB from a spe- 
cific location in the sky - GW radiometry using a network 
of detectors. We find that the angular resolution essen- 
tially depends on the effective GW bandwidth and the 
linear size of the network. In this paper we will restrict 
ourselves to the network of the two 4km LIGO detec- 
tors at Hanford (LHO) and Livingston (LLO). For the 
purposes of our analysis, we take their noise curves to 
be identical with the LIGO-I design power spectral den- 
sity [14]. Our future plan is to include VIRGO and other 
detectors around the world in the numerical implemen- 
tation of this analysis. 

The basic statistic is the cross-correlation between the 
data from a pair of detectors. In order to 'point' the 
pair of detectors at different locations one must suit- 
ably delay the signal by the amount it takes for the 
GW to travel to both detectors corresponding to the 
source direction. This delay will be a function of the 
source position and will vary as the Earth rotates. Using 
the delay allows the detectors to sample the same wave- 
front from the source. The cross-spectrum formulation 
has been carried out in [15, 16]. Methods for searching 
for isotropic backgrounds [17] using the cross-correlation 
and for anisotropies using spherical harmonic decomposi- 
tion [18] have been devised. Efforts have also been made 
to devise methods to measure the spherical harmonic mo- 
ments of the SGWB anisotropy using a network of ground 
or space based detectors [19, 20]. Here we focus on a spa- 
tially resolved search and the final goal is to make a map 
of the true SGWB sky. We achieve this goal by pixelizing 
the sky, that is, we use a pixel basis. 

The advantage of a spatially resolved search is seen im- 
mediately if we examine the so called overlap reduction 
factor, which partially determines the fractional power of 
source spectrum the search filter will receive at different 
frequency bands. The overlap reduction factor, normally 
denoted by j(f) in the literature [17] for the isotropic 
unpolarized background, becomes a time-dependent fac- 
tor j(fl, f,t) for the spatially resolved search. For the 
LIGO detectors, j(f) quickly reduces to zero beyond few 
tens of Hz, while j(Ct, /, t) has infinite bandwidth. So the 
bandwidth of the spatially resolved search is essentially 
detector bandwidth limited. This is typically valid for 



a network of detectors and therefore important from the 
point of view of the sensitivity regime of GW detectors 
which lies in this region. 

As in radio-interferometry, the correlation statistic so 
constructed produces a 'dirty' map where a point source 
does not produce a point image, but one that is smeared 
by a beam response function (beam, for brevity). The 
'cleaned' GW sky map is obtained from the measured 
cross-correlation statistic by deconvolving the beam. In 
other words, to obtain the GW power from each direction 
in the sky one needs to solve an integral equation where 
the measured power (data) is a convolution of the actual 
power with a kernel (beam) . In order to understand the 
structure of the beam we carry out a numerical and an 
analytical study using the stationary phase approxima- 
tion (SPA). We find that at low declinations (latitudes) 
of a point source, the kernel essentially has the shape of 
a figure '8' with a bright spot at the intersection. The 
bright spot is at the location of the point source. The fig- 
ure of '8' continuously changes and bifurcates into a 'tear 
drop' as the point source moves to higher declinations. 
The declination at which this bifurcation occurs is de- 
termined by the half-angle of the cone traced out by the 
vector joining the two detectors. For the LIGO detectors 
this declination is about 26° . The size of the bright spot 
or effective sky patch, defined by a certain percentage 
of reduction in the beam response function, say 50%, is 
determined by the inverse of the band-width divided by 
the light (GW) travel time between the detectors. Con- 
sidering a broadband source and LIGO detectors having 
kHz bandwidth with 10 ms light travel distance between 
them, the angular size is about 5° in radius. We find 
that these results agree very well with those obtained by 
applying singular value decomposition to the kernel ma- 
trix; the eigenvalues fall off steeply after a certain point 
which determines essentially the number of 'degrees of 
freedom' of the kernel matrix and thereby the size of the 
sky patch. 

We employ the maximum likelihood (ML) approach 
for deconvolving the sky map. The integral equation 
for a discrete pixelised sky leads to a set of linear al- 
gebraic equations. Several deconvolution algorithms ex- 
ist in literature for solving such a problem. However, 
because of the broad similarity of our present problem 
with the cosmic microwave background (CMB) analysis, 
we have opted for techniques that have been successfully 
applied to deconvolve CMB sky maps. Moreover, the 
ML approach provides a framework to study the SGWB 
anisotropy in other basis of interest, e.g., the spherical 
harmonic basis. We find the ML estimate by employing 
the conjugate gradient method. To verify our method, we 
apply the analysis on simulated sky maps mapped by a 
GW radiometer consisting of LIGO pair of detectors LHO 
and LLO. We generate 'realistic' colored Gaussian noise 
corresponding to LIGO-I design sensitivity curve [14] de- 
tector noise and embed various simulated sky maps of the 
GW stochastic background into the noise. We demon- 
strate that the true sky maps can be recovered satisfac- 
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torily. 

The paper is organized as follows: In section II we 
briefly review the GW radiometer concepts and obtain 
an expression for the GW radiometric cross-correlation 
signal, which is then optimized for the maximum signal- 
to-noise (SNR). In section III, we set up the integral 
equation that must be solved in order to obtain the true 
sky map from the data. A time-frequency analysis is 
performed and the directed optimal filter is derived for 
anisotropic searches. Moreover, a stationary phase anal- 
ysis is presented which provides us with the understand- 
ing of the kernel or the point spread function. In sec- 
tion IV, we describe the maximum likelihood approach 
and the conjugate gradient method and apply it to sim- 
ulated data to test the efficacy of this method. Later 
subsections outline the extension of the GW radiometer 
analysis to incorporate multiple baselines obtained with 
a network of detectors and the extension of the GW ra- 
diometer to search for polarized SGWB. The numerical 
implementation of the method is described in section V. 
We conclude in section VI. 




Detector 1 Detector 2 



FIG. 1: Geometry of an elementary radiometer. Above, 
Ax(i) is the separation or baseline vector between the two 
detectors; as the Earth rotates, its direction changes, but its 
magnitude remains fixed. The direction to the source ft is also 
fixed in the barycentric frame. The phase difference between 
signals arriving at two detector sites from the same direction 
is also shown. 



II. GW RADIOMETER EMPLOYING EARTH 
ROTATION APERTURE SYNTHESIS 

A. The principle of a radiometer 

Radiometry or aperture synthesis is a well known tech- 
nique in radio astronomy and CMB experiments. The 
idea is to point a pair of detectors separated by a baseline 
to a desired direction in the sky by introducing an appro- 
priate time-delay between their data-streams. This delay 
corresponds to the difference between the times of arrival 
of a GW signal if it were to arrive at those two detector 
sites from that direction. For a given source in the sky, 
this delay will change as the baseline orientation changes 
due to the rotation of the earth. The cross-correlation 
of the data from the two detectors, appropriately time- 
delayed, would cause potential GW signals arriving from 
the chosen direction to interfere constructively. Whereas 
signals from other directions will tend to cancel out be- 
cause of destructive interference. This principle of Earth 
rotation aperture synthesis, which is well-known in radio 
astronomy, could very well be used in GW astronomy 
using pairs of GW antennae. Figure 1 illustrates the 
principle on which the GW radiometer works. 

We consider the Celestial Equatorial frame whose ori- 
gin coincides with the centre of the Earth. The axes are 
defined as follows: For a fixed but arbitrarily chosen time 
t = 0, the x-axis is directed towards the intersection of 
the equator and the longitude <p = 0, which can be taken 
as the Greenwich meridian; the z-axis is directed towards 
the North Celestial Pole and the y-axis is chosen orthog- 
onal to the previous two axes forming a right-handed 
triad. The Earth rotates in this frame with the angular 
velocity lo e — 2n/(l sidereal day) ~ 7.3 x 10~ 5 radi- 
ans/sec oriented along the z-axis. The two detectors are 



at locations xj (where / = 1,2), and the baseline vec- 
tor joining the two sites is Ax := X2 — Xi. f2 is the 
unit vector in the direction to the source, which is fixed 
in the barycentric frame. The baseline Ax rotates with 
the Earth, but its magnitude |Ax|, which is the distance 
between the detectors, remains constant. The map of 
the SGWB can be constructed by performing this syn- 
thesis for each location in the sky, patch by patch. An 
approximate size of the patch or resolution of the GW 
radiometer and number of patches required to cover the 
full sky can be estimated from the following simple ar- 
gument. This simple model produces a fringe pattern 
with resolution given by the standard formula for the 
central width A8 ~ Agw/(| Ax| sin 6*) where Aqw is the 
GW wavelength and 9 is the angle of incidence. This 
is, however, a naive estimate. A better estimate of the 
resolution would follow from considerations involving the 
pixel-to-pixcl Fisher information matrix. In which case 
the solid angle resolution will scale inversely proportional 
to the square of the SNR. 

In this paper we will consider a GW radiometer made 
using the two 4 km LIGO detectors LHO and LLO in 
the United States. However, this analysis is equally well 
applicable to any other baseline or a network of baselines 
involving other detectors such as VIRGO, TAMA, GEO 
etc. The light travel time between the LHO and LLO is 
~ 10 ms, which at a bandwidth of 1 kHz gives a reso- 
lution A9 of few degrees. This resolution implies that a 
few thousand patches or pixels are required to completely 
cover the full sky. This estimate, in turn, will be useful 
for assessing the computational cost and numerical com- 
plications in handling large matrices for obtaining the full 
skymap. We make these estimates more robust below. 
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Basic Framework and the Statistical Properties 
of the Data 



We first set up the notation and framework needed for 
investigating the problem. 

In the transverse traceless (TT) gauge the metric per- 
turbation hij(t,x) in the equatorial frame can be ex- 
panded in terms of plane waves of the two polarizations 
+, x in the following way: 



fcy(i,x) 

where a 
form; the 

Here A = 
tion over 
spherical 
given by 



df / dn/u(/,n)eA(n)e 



2iri/(t+f2-x/c) 



(2.1) 

tilde on a variable denotes its Fourier trans- 
complex Fourier amplitudes satisfy the relation 
= hA.(—f,£l) owing to the reality of hij(t,x). 
-- {+, x} is the polarization index and summa- 
the repeated index A is implied. In terms of the 
polar coordinates (6, cf>), the source direction is 



ft = cos 6 sin ( 



sin sin f 



COS t 



(2.2) 



and, hence, the infinitesimal solid angle along the direc- 
tion ft is dft = sin 9 &Q d<f>. Note that the wave propaga- 
tion direction is — ft. The polarization tensors e (ft) are 
defined by the following expressions: 



e + («) 



eg ® e e - 
eg <g> e<* + e„ 



1 , 
) e 9 , 



(2.3a) 
(2.3b) 



where 



eg = cos <j) cos # x + sin <f> cos # y — sin 9 z , 
= — sin x + cos <f> y , 

are the two orthonormal basis vectors on a two-sphere. 

The statistics of the GW signal can be best described 
in the Fourier domain: If we assume the signal to be 
stochastic and uncorrelated in the two polarizations [50], 
different frequencies and different directions, then the 
Fourier components of the GW strain obey, 



(h* A (f,n) M/'.n')) 
= 8 AA ,s(f-f)s 2 (n 



tl')V A {tl)H{f), (2.5) 



where Va(&) is proportional to the strength of the 
SGWB in the direction ft and H(f) is the two sided 
GW source Power Spectral Density (PSD). The inter- 
pretation of the quantity Va^) can be made apparent 
by relating it to the specific intensity or brightness [21] of 
GW /gw(/, ft)- Specific intensity is defined as (c times) 
the incident energy density per unit frequency interval 
per unit solid angle, or, equivalently, specific intensity is 
the (normally incident) flux per unit solid angle. Fol- 
lowing the convention commonly used in SGWB anal- 
ysis, if the incident GW energy density is expressed in 
the units of critical energy density for a flat universe 



p c := 3c 2 H 2 / (8irG) , where Hq is the Hubble constant 
at the current epoch and G is the Newton's gravitational 
constant, it can be shown that [22]: 



W/>«) 



47T 2 C 

3H 2 



f 2 H{f) T+(Ct) + p x (fi) . (2.6) 



Therefore, Va(£1) can be interpreted as the specific in- 
tensity of SGWB for the corresponding polarization (up 
to a certain proportionality constant). 

In general, we cannot separate H(f) from Va(&), be- 
cause the frequency power spectrum H(f) could also 
depend on the direction ft. A more general treatment 
would use a quantity like V A (£l, f) which describes both 
frequency and angular distribution of SGWB power to- 
gether. But for a small enough bandwidth the assump- 
tion may be justified, as the signal is expected to have a 
smooth profile of the power spectra. Further, it should be 
noted that V a is actually a second rank tensor and should 
be represented by two indices as Vaa', but because we 
assume that the two polarizations are uncorrelated the 
quantity V+ x is zero and thus absent. To avoid unnec- 
essary indices and with a slight abuse of notation, we 
therefore write V a with a single index instead of V aa' ■ 

We consider two GW detectors located at x/ (t),I = 
1,2. The detector frames are denoted by the coordinates 
(Xj,Yi, Zj) where the arms of the respective detectors 
lie along the (Xi,Yj) axes. Then the detector tensors di 
are given by: 



'Y, 



(2.7) 



The factor | is inherited from the geodesic deviation 
equation. Owing to the Earth's rotation, the detector 
coordinates and di are functions of t. Thus in matrix 
form: 



di(t) = Tl T {t) ■d I (t = 0) -n(t), 



(2.8) 



where TZ(t) is the rotation matrix describing a rotation of 
angle LUEt around the Earth's rotation axis, namely, the 
z-axis. Here we is the Earth's sidereal angular velocity 
w 7.3 x 1CT 5 radians/sec. 

The strain hj(t) in the I th detector is given by: 



h I (t) = h lJ (t,x I (t))di j (t). 
We define the antenna pattern functions as 
Ff(t;fl)=ef j (n)d i I j (t) 



(2.9) 



(2.10) 



for a wave incident from the direction f2. 

Contracting Eq. (2.1) with the detector tensors dj, 
the signal amplitudes can be expressed in terms of the 
antenna pattern functions as: 



hj(t) 



df / dClh A {f,Cl)Fj A {t;Vl)e 27Ttf{t+(l -* /c) . 
Js 2 

(2.11) 
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The baseline Ax at any time t (in matrix form) is given 
by: 

Ax(t)=K T .Ax(i = 0). (2.12) 

The detector output sj(t) is a sum of the GW signal 
hj{t) and noise nj(t) 

si{t) = h I {t) + m{t). (2.13) 

Statistically, the gravitational- wave strains hi(t) are un- 
corrected with the detector noise n.j(t); that is, the fol- 
lowing four correlations, in the time domain are zero: 

(hi(t)nj(lf)) -0 J, J =1,2, (2.14) 

where t, t' are any two time instants. We also assume 
that the noise in different detectors is uncorrelated; that 
is, (ni(t) ri2(t')) = 0. This assumption is not unreason- 
able when the detectors are widely separated. Thus the 
only possible correlation is between h\(t) and /^(i')- The 
statistic that we construct in the next subsection is based 
on this fact. 



C. The Cross-correlation Statistic for the Directed 
Search 

Stochastic GW signals inherently can arrive from any 
direction with any amplitude. Moreover, they are char- 
acterized by the statistical expectation values of energy 
density. The noise in the two detectors are assumed to 
be essentially independent. In this situation the cross- 
correlation of the data from the two detectors is an ap- 
propriate statistic for detecting and observing stochastic 
GW. In order to optimize the SNR, the cross-correlation 
statistic for the directed search involves a direction- 
dependent filter function Q(t,Q.;t' , t"), which connects 
sidereal time t' of one detector's data to t" of the other 
detector's data to match the phases of the GW strains in 
the detectors. The filter does a inverse noise weighting 
using the PSDs of the two detectors in order to suppress 
noisy frequency bands and enhance the SNR by assign- 
ing relatively large weight to the sensitive regimes of the 
detectors and the bands where the source is expected to 
emit more power. (In general, Q will depend on Va{&)- 
However, for the directed search that we envisage here, 
as we will see later, the Va{&) are delta functions and, 
therefore, the filter Q only depends on tt.) 

The sampling interval with which the data are sampled 
is determined by the Nyquist frequency of the stochastic 
signals of interest, and can be well below a millisecond, 
corresponding to several kilohertz. While it is possible 
to compute the filter function on a time segment of this 
sampling interval, it is erroneous to do so for the physical 
problem at hand. This is because the signals at the two 
detectors will be incoherent on time-scales smaller than 
the light-travel time, td, of that baseline, which is at most 
a few tens of milli-seconds for ground-based detectors. 



Thus td sets the lower limit on the duration of the time 
segment on which the filter function should be computed. 
The upper limit, r, is set by the smaller of the timescales, 
on which the data are stationary and the timescale on 
which the detector orientation changes appreciably. We 
thus divide the data into time segments, At, such that 
td <C At <C t. The time segments used currently in 
LSC data analysis vary from 32 to 192 seconds, and are 
consistent with these limits. 

The final statistic for the full observation time T is ob- 
tained by linearly combining the cross-correlations over 
the smaller time-intervals as a weighted sum. The filter 
is optimized for each time segment, I k = [tfe — At/2, t k + 
At/2], at sidereal time t k and the statistic for the k th 
segment is given by: 

AS(t = t k7 Cl) = AS k (Cl) 
= [ dt' [ dt" Sl {t')s 2 {t")Q(t k , tl;t',t"). (2.15) 

The final cross-correlation statistic S for all the n = 
T/At sidereal time bins can then be obtained by com- 
bining the AS k as a weighted sum as follows: 

n 

= AS fc) (2.16) 

fc=i 

where, the n quantities w k (Ct) are to be chosen so that 
the SNR for the statistic S(ti) is maximized. We denote 
the SNR of S by ps = ps/o~s, where ps and 175 are the 
mean and standard deviation of S respectively. We use 
normalized weights J2k=i w fe = 1- 

The Us, <Js and the SNR pg are given in terms of the 
means pk = (ASk) and variances a k = (AS|) — (ASk) 2 
of the individual mutually uncorrelated time segments as 
follows: 

n 

MS = ^W k p, k , 
fc=l 
n 

2 \ ^ 2 2 

a s = 2^ w k a k^ 

k=l 





n 




n 


PS = 


^2w k p k 


/ 


\ ^ 2 2 

}^ w k Ok 




_fc=i 




.fe=i 



The maximization is achieved elegantly by invoking the 
Schwarz inequality. We consider an n-dimensional Eu- 
clidean space equipped with the usual scalar product 
in which w is a vector with components w k . We de- 
fine k and p having the components n k = w k a k and 
Pk = Pk/&k respectively. Then, we may write ps = & ■ p 
where the vector k is a unit vector in the direction of 
k. pg is maximized when k points in the direction of 
p. Thus, w k cx p k /a k , the proportionality factor being 
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Efe=i Mfc a k 2 ] • W e have the results: 



S = 



En -2 
fc=i Mfc CT fe 



AS 



= Ip1= 



(2.18) 



(2.19) 



Further, we are free to choose the normalization of the 
filter Q. We choose the normalization such that all fi^ 
are made identical, equal to fi, say. This simplifies the 
statistic S and its SNR p$- 



S = 



v^'i -2 




(2.20) 
(2.21) 



which are their final forms. Normally the GW signal is 
expected to be weak so that even after its integration 
over At it is still smaller than the noise, that is, <C 
CTfe. Therefore the variance in each segment obeys a\ ~ 
(AS 2 ). The signal however builds up when we integrate 
over all the time segments during the observation time. 
Secondly, from Eq. (2.20), it is evident that noisy time 
intervals contribute less to the statistic S leading to its 
optimal character. 



III. THE CONVOLUTION EQUATION FOR 
THE SKY MAP 



Since the GW power spectra and the detector noise are 
modeled in the frequency domain, it is convenient to for- 



mulate the whole analysis in that domain. However, since 
the detector output is a time-series, it is pertinent to ask 
over what time-duration must one compute their Fourier 
transforms. As discussed in the last section, typical ac- 
ceptable segment sizes are a few tens to a few hundreds of 
seconds. Fourier transforms computed over such "small" 
(vis a vis the total observation) time scales are termed 
as short-term Fourier transforms (SFTs). Each SFT be- 
comes a function of time t as well because t is essen- 
tially the identifier of the segment. Thus we have a time- 
frequency representation of the data. As we shall see, 
this representation is most suitable for further analysis 
and has also been used in previous literature [17, 18, 19]. 



Time-Frequency Analysis of the Signal and 
Noise 



The (approximate) SFT of a segment of detector out- 
put can be defined as [18], 



i-t+At/2 

si(t;f) := / dt's^e- 2 ^' . (3.1) 

Jt-At/2 

We retain here the convention of using 'tilde' over a sym- 
bol to denote Fourier transform - the distinction should 
be evident from the context. Most importantly, by taking 
the inverse Fourier transform. 



/ 

J — < 



dfsi(t; f)e 



iizift 



rt+At/2 

df e 2 " /4 / At' si (f) e~ 27rift ' = sj (t) 

i Jt-At/2 



(3.2) 



the exact time series segment sj (t) can be recovered. The 
same notation will be used for several other quantities as 
well in this analysis. 

While constructing the statistic, it was assumed that 
the true GW strains hj(t) in the detectors are correlated, 
but that the noise streams are uncorrelated. The expres- 
sion for the correlation between the GW strains in two 
detectors will be derived here, which is necessary for the 



derivation of the optimal filter in the next subsection. 
The detector parameters can be approximated to be sta- 
tionary over the period of the time segment. Hence we 
may consider the quantities dj and x/ to be nearly con- 
stant over a time segment and regard them as functions 
of the time t labeling the time segment. The SFT of the 
GW strain in detector / over a time segment is given by: 
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hi(t;f) 



roc r ~ 

dt' / df / dm A {f',n)Ff(ci,t)e 

t-4* J~oo Js 2 



27Ti f't'-ft'+f 



s 2 



dOF/(n,t) / d/ 



(3.3) 



r 



where summation over A is implied. The <5a*(/) is the 
finite time delta function (sine function) defined by, 



-At/2 "7 



(3.4) 



The finite time delta function <5a*(/) behaves as the Dirac 
delta function in the limit At — > oo, but has the 
property 5a* (0) = At. Hence for a large time segment 
At the SFT from detector / takes the simple form: 



hi(tj) 



dnh A (f,Cl)Ff(Cl,t) e 



27ri/n-xj(t)/c 



(3.5) 



The important result of this subsection is the expec- 
tation of cross-correlation between the SFTs of time seg- 
ments of detector outputs at time t from the two detec- 
tors 1 and 2, which can be obtained from Eq. (3.3) and 
Eq. (2.5) as: 



(hl(t, /) h 2 (t, /')) = e 2 ^f~^ / df H(f") 7 (f, /"; AQ, dV A )5 At (f - /) 6 At (f" - /') 



7(t,/;Afi,dP A ) 



dn 



AQ=S 2 



F+(n 1 t)F+(n,t)r+(n) + F?(si,t)F?(si,t)p x (si) 



_,27ri/n-Ax(t)/c 



(3.6) 



(3.7) 



r 



The general overlap reduction function defined by 
Eq. (3.7) is a generalization of the usual overlap reduction 
function for the isotropic SGWB case, first constructed 
by Nelson Christensen [15] and formally written in a 
closed form by Flanagan [16]. In this case it is a more 
complex object. Besides the frequency, it is also a func- 
tion of the segment time t. It also depends on V A (&) 
which is integrated over the full sky S 2 . Thus it is a 
functional of V A (£l). In general, when we construct our 
directed filters, the integral will be restricted to a small 
patch of the sky, Af2 c S 2 . (In principle, Af2 can be any 
(measurable) subset of S 2 .) Accordingly, V A (tl) plays 
the part of a weight function over the sky and we may 
define the measures: 



dV A = V A (Cl)dCl . 



(3.8) 



Thus, the 7 in general becomes a functional of the SGWB 
power in both polarizations coming from the patch Afl. 
We therefore separate the function arguments, t and / 
from the non-function arguments, Ail and V A , by a semi- 
colon. Finally, the exponential term in Eq. (3.6) be- 
fore the integral is just the time-shift term in the Fourier 
transform of the segment at time t. 



In the limit of a large time segment, Eq. (3.6) takes 
the simple form: 

(ht(tj)h 2 (tj')) = 5(f-f)H(f) 1 (t,f-S 2 ,dV A ). 

The advantage of expressing (h*(t,f)h2(t,f)) by 
Eq. (3.6) can be readily realized if we put / = /'. In this 
case, the correlation given by Eq. (3.9) diverges in the 
limit At — > 00. But, in practice, At is finite, and hence, 
we expect a finite value for the correlation. Eq. (3.6) 
lets us compute that finite value of (h\(t, f) h 2 (t, /)) at 
/ = /' . We use the large At limit and replace one of the 
finite time delta functions 5/\ t (f" — f) m the integrand 
of Eq. (3.6) by the Dirac delta function S(f" — /), while 
treating the other <5a*(/" — /) as a normal function and 
put <5 At (0) = At. We get, 

(hl(t,f)h 2 (tj)) = AtH(f) 1 (tJ;S 2 ,dV A ). (3.10) 

This formula was derived by following the same pro- 
cedure as prescribed in [17], so, not surprisingly, for 
isotropic backgrounds our result matches the formula ob- 
tained in [17]. This result is important for injecting test 
signals in the detector output [23]. 
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Next we describe the properties of detector noise in a 
finite time segment. 

The noise in the segment labeled by t is a time series 
ni(t) in a detector /. Its SFT is given by: 

pt+At/2 

nj(t;f) := / dt n,(t) e" 2 ^*. (3.11) 

Jt-At/2 

We take the mean to be zero: (ni(t)} = (n/(i;/)) — 0. 
Since rij(i) is real, its SFT obeys the reality condition, 
n*j(t; f) = ni(t;—f). The noise in a detector is un- 
corrected with the noise in another detector and with 
the GW signal, i.e., (m(t) n 2 (t')) = (fci(t) n 2 (t')) = 
{n l {t)h 2 {t')) = [see Eq. (2.14)]. These relations also 



In the limit of large length of the time segment, we arrive 
at the usual formula 

(nK*;/)nz(t;/0) = \s{f - /') Pz(t; |/|). (3.14) 

Again, the advantage of using Eq. (3.13) in expressing 
(nj(t; f)ni(t; /')) becomes evident when we set / = /'. 
The usual formula, Eq. (3.14), involving Dirac delta func- 
tion diverges, whereas, in practice, that expression is ac- 
tually finite because of a finite observation time. How- 
ever, in Eq. (3.13), if we replace one finite time delta 
function by the Dirac delta function and treat the other 
as a normal function we obtain a finite result, 

(\nj(t;f)\ 2 ) = ^Am(t;|/|). (3.15) 

This is used in our work when generating simulated noise 
for our test analyses. 



B. Optimal Filters for Anisotropic Searches and 
the Directed Search 

The aim of this subsection is to construct an optimal 
filter to maximize the SNR of the cross-correlation statis- 
tic over the small time-segments. We essentially general- 
ize the analysis presented in [17] for isotropic background 
search to the case of anisotropic background search. 

The optimal filter depends on the theoretical model of 
the SGWB, that is, on P A (£i) and H(f). First we will de- 
rive the filter for the general case of the anisotropic search 
and then specialize to the directed search. Our first goal 
is to compute the SNR p{t) of the statistic AS(t) over 
the time segments of length At at time t. 



hold for their corresponding SFTs. The length of the 
time segment is usually kept a few tens of seconds long, 
over which the detector noise can be regarded as station- 
ary. Thus (nj(t') ni(t")} is a function of t" — t' , provided 
both t',t" are in the same segment centered at time t. 
Then, using the fact that ni(t) is real, we have, 

1 r°° ,, , 

(nz(t')nz(t")) = ~]_ d/^ftl/Oe 2 **** "°. 

. (3 ' 12) 

where Pi(t; /) is the one-sided noise PSD. This noise 
PSD is also a function of time t as detector noise is non- 
stationary. The correlation between the corresponding 
SFTs can be easily obtained from the above relations: 



(3.13) 
I 

The equation for the general case is the generalization 
ofEq. (2.15): 

AS(t)= [ dt' [ dt" Sl (t')s 2 (t")Q(t;t',t"), 
Ji(t) Ji(t) 

(3.16) 

where I(t) is the interval [t - At/2,t + At/2}. Here we 
have suppressed the model dependence of Q. Assum- 
ing that the noise in both detectors and the earth are 
stationary within the duration of each time segment, we 
may write Q(t;l/,t") — Q(t;t' — t"), which allows us to 
expand the filter in terms of its SFT Q(t, f) as, 

/oo 
Afe ^W-t") Q{t j), ( 3.i 7) 
-00 

Thus in terms of SFTs the statistic can be expressed as, 

/oo 
dfsl(t;f)s 2 (t;f)Q(t,f). (3.18) 
-OO 

The mean of the statistic AS(t) is: 
H(t) := (AS(t)) = ^ df(Sl(t;f)s 2 (t;f))Q(tJ). 

J — OO 

_ (3.19) 

Replacing si by hi within the ensemble average, since the 
relevant correlations of signal and noise are zero except 
for (hl(t; f)h 2 (t; /)), we obtain from Eq. (3.10), 

/OO 
dfH(f) 1 (tJ;S 2 ,dV A )Q(t,f). (3.20) 
-00 

Here the H and V a are the actual quantities pertaining 
to the SGWB source. The corresponding quantities of 
the theoretical model are hidden inside the filter Q. 



(n}(t; f) n,(t; /')) = \ f d/"P 7 (t; |/"|)M/" - /) M/" - /') • 
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The variance after a routine but fairly involved calcu- 
lation is obtained as: 

a 2 (t) = ([AS(t) - (AS(t))} 2 ) 



At 

T 



d/P 1 (t;|/|)P 2 (t;|/|)|Q(i,/)| 2 .(3.21) 



From these results the SNR p(t) for the segment at time 
t can be computed. However, we need to maximize the 
SNR over each time segment. Here again we follow the 
prescription presented in [17] - we invoke the Schwarz 
inequality. To this end it is convenient to define a scalar 
product of two functions A and B on each time segment, 
labeled by t, as, 



(A,B)(t) 



dfP 1 (t;\f\)P 2 (t;\f\)A*(t;f)B(t;f) 



(3.22) 

The norm of a function A(t) is defined as || A || 2 = (A, A). 
Then the mean and variance are given in terms of the 
scalar product as follows: 

^ ) = Ai fm^ ig(i) W 



a 2 (t) 



V Pxfa |/|) P a (t; |/|) 
-At || Q || 2 (t) . 



(3.24) 



The SNR pit) is just the ratio p(t)/a{t). The SNR is 
maximized when the "signal" and the "filter" vectors are 
parallel, which happens when, 



Q(t,f) = X(t) 



g(/) 7 *(t,/;5 2 ,dP A ) 
P(i; l/l) P 2 (t; 1/1) ' 



(3.25) 



where X(t) is a (real) proportionality constant for the 
time segment. It is a function of the segment time t. 
This function will be chosen so that the SNR of S for the 
full observation is maximized. 

For the optimal filter, the expression for the mean 
given by Eq. (3.20) simplifies to: 



n(t) 



At II Q II 2 (*) 
X(t) ■ 



(3.26) 



We exploit the freedom of choosing X(t) by setting p{t) 
1 for each time segment. This immediately gives: 



X(t) = At||Q|| 2 (t), 



a\t) = 



\{t) 



(3.27) 
(3.28) 



We further require to find X{t) explicitly. This is done by 
computing || Q || 2 . We find || Q \\ 2 — X 2 P^ W where, 

M 2 (f)h(tJ;S 2 ,dV A )\ 2 



Pnw(1) 



df- 



P 1 (t;\f\)P 2 (t;\f\) 



(3.29) 



The integral on the right-hand side of the above equa- 
tion is positive definite: Its integrand contains the fourth 



power of the GW amplitude. Therefore, it may be de- 
noted by the square of a real quantity Pnw, which is 
determined by the GW power accessible to the network 
of the two detectors. The above equation gives the nor- 
malizations: 



Q(t) || 2 - |A<P w (t)]- 2 , 



(3.30) 



The optimal statistic is then easily obtained as in 
Eq. (2.20). We replace the sum in that equation by an 
integral over the segment time t: 



S = 



J AS(t)<j- 2 {t)dt 



(3.31) 



/<T- 2 (t)dt ' 

where the integration is over the observation time (which 
could consist of disconnected time intervals). 

The filter given by Eq. (3.25) is the general optimal 
filter to search for any anisotropic SGWB, which, unfor- 
tunately, requires the knowledge of H(f) and V A (Cl). In 
practice, we do not have an exact a priori model for H(f) 
and V A (£l), which anyway we are trying to measure! We 
then must use models for those quantities, H'{f) and 
V' A {tl), to search for different anisotropic backgrounds. 
These models will be used to construct the model- 
dependent overlap reduction function j(t, /; Af2, dV' A ) 
and the filter Q. 

The angular power distribution for only one unit point 
source on the sky in the direction f2 with equal power in 
both the polarizations can be expressed as: 

V A (£l') = <*(«'- ft). (3.32) 

This immediately simplifies the expression for the overlap 
reduction function because now the integral in Eq. (3.7) 
simplifies owing to the delta-function, and 7 becomes a 
function of f2 as well. Therefore, we have: 



7(n,t,/) = T(n,t)e 2mf{1 - Ax(t)/c 



r(n,t) 



A 



(3.33) 
(3.34) 



Unlike the time-independent overlap reduction function 
of the isotropic SGWB case, the direction-dependent 
overlap reduction function, j(fl, t, /), accepts power from 
all the frequencies and in fact has infinite bandwidth in 
the limit of vanishing pixel area. So the bandwidth of the 
filter Q in this case would only be limited by the band- 
width of the detectors through the coefficients Pj (t; /). 

If, instead, the source has a finite spatial extent, the 
bandwidth would be limited, because the integral in 
Eq. (3.7) would have to be performed over the solid angle 
Afi subtended by the source. If one takes a small patch 
of the sky of size (A9, A<f>) around some fixed (source) 
direction f2 := (0O)<fo)> it is easy to show that 7 is es- 
sentially a product of sine functions in (A8, A<f>). In fact: 



a£i e 27ri/(n-n ).Ax/c 

. nfAdAx.ee . 7r/A0Ax.e 
|AS2| sine sine ,(3.35) 
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where |Afi| = sm9A8A(f> is the solid angle subtended by 
the patch Ail. In the integral, the factor T(tt, t) remains 
nearly constant. The sine functions go to zero when their 
arguments reach 7r radians. Taking this definition as the 
bandwidth and taking |Ax|/c ~ 10ms for the two LIGO 
detectors, the bandwidth is about 750 Hz for a square 
source of side 10° on the sky. 

If there were known models for the anisotropic SGWB, 
the optimal filter for the general anisotropic case would 
have included V' A (Cl) and we would perform a full sky 
search for an anisotropic background. However, no rea- 
sonable model for the anisotropic SGWB sky exists in lit- 
erature and, so, blind estimations are currently the only 
possible alternatives. 

Directed search is one blind estimation approach, 
where the strength of each point (pixel) of the sky is 'ob- 
served' using a direction dependent filter, assuming that 
the other points on the sky do not contribute towards 
the observed value. So, in the directed case Q becomes 
a function of f2: 

where X(tl,t) is the normalization constant, which now 
varies from pixel to pixel. It is given by Eq. (3.30), but 
with the 7 in the expression for P% w in Eq. (3.29) re- 
placed by the simpler T of (Eq. (3.34)). The directed 
filter given in Eq. (3.36) is an optimal filter if there is 
a single point source in the direction f2 and no sources 
elsewhere in the sky. If there are other sources in the sky, 
as in a general anisotropic background, the filter becomes 
suboptimal as it stands. However, we use the above filter 
to make a "dirty" map of the sky, which is a convolu- 
tion of the the actual anisotropic background with the 



beam function and contains additive noise. We intend 
to extract information about the true background by the 
process of deconvolution. The convolution equation will 
be derived in the next subsection. 

The working principle of the above filter is evidently 
similar to the earth rotation aperture synthesis often 
used in CMB and radio astronomy to make map of a cer- 
tain portion or the whole sky. The phase lag between two 
detectors, separated by a distance Ax(i), in receiving a 
plane wavefront from a certain direction f2, as shown in 
Figure 1 , is compensated in the filter via the phase factor 
exp[2Triftt ■ Ax(i)/c]. As the earth rotates this factor 
adjusts, so that, waves from the given direction are 
coherently added, while the waves from other directions 
tend to cancel out. Note that, we did not introduce the 
phase factor by hand, it appeared automatically through 
the process of the maximization of the SNR. Though 
the whole radiometer analysis is based on this principle, 
the idea is clearly realized in the directed search analysis. 

C. The Integral Equation for the Directed Search 

In this subsection we set up the convolution equation, 
which is an integral equation for the statistic S(tt) - the 
dirty map of directed search. The kernel of the integral 
equation consists of beam functions that we define below. 
The goal is to obtain the power in both polarizations 
Va^), given the statistic S(tt). To this end we take 
the expectation value of the statistic S(tt) in Eq. (3.31) 
and use Eq. (3.20) for the expectation value of AS* inside 
the integral sign. We must also use expressions for the 
source overlap reduction function from Eq. (3.7) and the 
directed filter from Eq. (3.36). The final form of the 
(noiseless) convolution equation is given by: 



s (Ci) ee (s(Ci)) = [ dfi' \B + (n,n')v+(n') + b x (Ci,Ci , )v x (Ci / )] , 

Js 2 L J 



(3.37) 



where the beam functions B (i~2,f2') are defined as: 



B A (n,Cl') = A(«) / dt 



Pi(*; |/|) iVt; |/|)' 



A _1 (n) 



1 

At 



\~ 1 (t,n)dt. 



(3.38) 
(3.39) 



The function H'(f) is the model power spectrum of the 
source we insert in the kernel and Af2 = tt — tt' . We 
measure (S(tt)} and from the kernels B (tt, tt'), we pro- 
pose to solve the integral equation for VA(tt). Physically, 



we may expect the power in both polarizations to be the 
same, that is, "P + (f2) = "P x (f2) = ^(tt), say, the kernel 
then is just the sum of the two individual kernels of each 
polarization, 
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B(n,n') = b + (Ci,Ci') + b x (Ci,Ci') 

= m I dt £ if mwrnm r( "' " r( "' " e " 2 " ,A " ■ (3A0} 



Our numerical deconvolution strategy is described in the 
next section. But before we do that we examine the 
kernel in Eq. (3.40) and try to understand it from a 
physical point of view. This will afford us some insight 
into the beam patterns associated with directed filters. 

It is worth noting that the beam function B(tt,tt') 
is not symmetric only due to the leading normalization 
factor A(i~2), which comes from the normalization of the 
statistic S(tt). We make use of this observation to intro- 
duce a symmetric kernel in section V, which is advanta- 
geous for numerical deconvolution. 

D. The Stationary Phase Approximation of the 
Kernel and its Singular Value Decomposition 

The GW radiometer beams are not pointed but have 
a spread out profile, which varies with sky position. 
Thus in order to make progress towards deconvolving 
the GW sky map we try to understand the beam pat- 
tern. We find that the Stationary Phase Approxima- 
tion (SPA) of B(n,ft), given in Eq. (3.40), yields 
useful results. It is essentially the exponential term 
cxp[— 2irifAtt ■ Ax(f)/c] containing the phase that de- 
termines the integral - the integrand constructively con- 
tributes when the phase in the integral in Eq. (3.40) 
is stationary. We also use the fact that rest of the func- 
tions in the integrand vary slowly with time, so that they 
effectively behave as constants as far as the integral is 
concerned. 

We obtain the beam function for a unit point source at 
tl = (6> , (f> ). We write Ail := Cl - Cl . Note that Af2 
may not be necessarily small; the points fi, tl Q can lie 
anywhere on the unit sphere. By performing a numerical 
computation for an observation time of one sidereal day, 
we find that for low declinations, the beam is shaped like 
the figure of "8", as shown in Figure 2(a), while as one 
goes higher in declination, the "8" smoothly turns into a 
"tear drop" . 

With the application of the SPA we can explain the 
shape of the beam. The integrand in the kernel usually 
oscillates rapidly, because of the exponential phase term, 
except when the phase is stationary. The kernel is a dou- 
ble integral over / and t and therefore the SPA must be 
carried out in two dimensions. Setting the first derivative 
of the phase with respect to both variables / and t equal 
to zero, we obtain: 

Aft • Ax(<) = 0, (3.41) 
An-Ax(f) = 0, (3.42) 



r 




-0.4 -0.2 0.2 0.4 0.6 0.8 1 

(a) Numerical beam pattern. 




(b) Theoretical beam pattern obtained by stationary 
phase approximation (SPA) 



FIG. 2: Illustration of the agreement between numerical 
and theoretical GW radiometer beam patterns at declina- 
tion + 12° for the LIGO detectors at Livingston and Han- 
ford (with white noise, upper cut-off frequency of 1024 Hz, 
_H"(/)=constant and observation time of one sidereal day). 
Same colormap has been used for both panels. Note that, 
the shape of the beam does not depend on the right ascension 
of the pointing direction for observation time of full sidereal 
day(s). 

where Ax(i) := dAx(i)/di. The detector separation 
vector Ax(<) rotates about the earth's rotation axis (z- 
axis in our coordinate system) with the angular velocity 
we- Geometrically Ax(i) traces out a right circular cone 
with z-axis as its symmetry axis [see Figure 3]. It is 
explicitly given by: 

Ax(t) = Ai?(sin0coswEijSin0sina>£i,cos0) 1 (3.43) 

where AR = |Ax(i)| is the constant distance between 
the detectors. As w^i ranges from to 2ir, Ax(t) traces 
out a cone with half angle O. The half- angle O of the 
cone, < < 7r, is given by, 

^ cos 9i — cos 8 2 

COS0 = =, 
y2[l — cos 6*i cos 02 — sin 0\ sin 02 cos(0i — (fo)] 

(3.44) 
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FIG. 3: The baseline formed by two detectors, Ax(t), traces 
out a cone as the earth rotates. A schematic diagram is shown 
here. The vector n cone (t) is normal to the baseline as well as 
the cone. 



where (0j,4>i) are the detector coordinates. For the 
LIGO pair of detectors, 6-27°. 

From Eqs. (3.41) and (3.42) it is clear that the phase is 
stationary when Ax(t), Ax(t) and Af2 form an orthog- 
onal triad. Since the unit vector normal to the baseline 
at any given time as well as to the cone 



Ax(t) x Ax(i) 



|Ax(t) x Ax(*)| ' 
the SPA condition is satisfied when: 



An = Afin, 



(3.45) 



(3.46) 



where Afl = |An| can take both positive or negative 
values. Since both n and tt(t) = CIq + Afi(f) are con- 
strained to lie on the unit sphere and thus both have unit 
norm, it follows from Eq.(3.46) that, the SPA solution 
tt(t) is a curve on the unit sphere given by [22], 



n(t) 



2[Ct 



com: 

(*)■ ( 3 - 47 ) 



The trajectory has been parameterized in terms of the 
sidereal time t. One can even obtain an approximate an- 
alytical expression for the beam function along the SPA 
trajectory using standard SPA techniques as [22]: 



B(fl(t),fl ) = A(ft(t))T(ft(t);t)T(ft ;t) x 



\ [z • Ax(t)] [z • («(t) - Ho)] 



.(3.48) 



As t is varied over a full sidereal day, the shaded figure 
of "8" is generated through Eq. (3.47) and Eq. (3.48) 
as shown in Figure 2(b). Clearly, SPA results match 
very well with the numerical beam pattern shown in Fig- 
ure 2(a). 

The case where Eq. (3.48) does not apply (though the 
analysis still remains valid) is when the detectors are at 
the same latitude, as the normal to the baseline cone, 
n cone (i), always remains parallel to the z axis, causing 
the denominator of Eq. (3.48) to vanish. In this case 
the whole SPA trajectory shrinks to a point, which is 
the image of the pointing direction about the equatorial 



plane, ft = tto — 2[tto ■ z]z. 
function at the image point, 



The value of the beam 



dtr(n ,t)r(n - 2[n -z]z,f) 



dt[T(n ,t)f 



(3.49) 

is also quite easy to compute. Therefore, a skymap pro- 
duced by such a baseline will be a superposition of the 
(blurred) true sky and its (differently blurred) reflection 
about the equatorial plane. In practice, a situation like 
this arises for the LHO-Virgo pair, as their latitudes are 
quite close, 46°27"N and 43°37"N respectively. 

The SPA solution given in Eq. (3.48) does not remain 
finite very close to the pointing direction, as the denom- 
inator vanishes. However, close to the pointing direction 
a better approximation is obtained by expanding the the 
phase term exp[— 2irif Af2 • Ax(f ) /c] up to the second or- 
der, which can also be used for more accurate modeling 
of the core of the beam. 

The stationary phase analysis also indicates an approx- 
imate resolving power of the radiometer. Since near the 
maximum of the beam function, the phase must not vary 
too much over the bandwidth, say no more than a ra- 
dian, this implies that the resolving power is determined 
by | An | ~ A/ 1 Ax |, where the A corresponds to the band- 
width / = A/ of the detectors. For the LIGO detectors, 
| Ax| ~ 3000km. If we take the bandwidth to be ~ 1kHz, 
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FIG. 4: Contour plot of GW radiometer beam pattern at dec- 
lination + 12° revealing the approximate size and the (nega- 
tive) side lobes of the beam for the LIGO detectors with white 
noise, upper cut-off frequency f u — 1024 Hz, (/)=constant 
and observation time of one sidereal day. The contours are 
drawn starting from the highest level of 0.9 with a difference 
in levels of 0.1. The beam falls by 1/e in between the 5th 
and the 6th contour (from the highest value), which is a good 
indicator for the beam size. Clearly, the narrowest beam size 
(along the "minor axis") is in reasonable agreement with the 
theoretical prediction of 6° ~ 0.1 rad. The beam becomes 
broader (not shown in the figure) for real detector noise and 
negative source spectral index, e.g., H(f) oc /~ 3 . 
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then A ~ 300km, so the radiometer resolution is ~ 0.1 
radians, that is, ~ 6 degrees, which we find is consistent 
with the numerically obtained beam profiles. A contour 
plot of the core of the beam of the radiometer formed 
by the LIGO pairs (with white noise and upper cut-off 
frequency of 1024 Hz) is shown in Figure 4. The plot 
confirms that the beam size in this case is 6° ~ 0.1 rad. 



the detectors is not as good, hence the achievable angu- 
lar sensitivity of the radiometer becomes relatively poor. 
The plot of singular values for the same detector pair 
both having LIGO-I goal noise PSD with an upper cut- 
off frequency of f u — 512 Hz is overlaid (dashed line) 
on Figure 5. Clearly, the number of degrees of freedom, 
which represents the amount of information content in a 
map, is less in this case. 
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FIG. 5: The solid line is a plot of the singular values of the 
kernel for the LIGO pair of detectors with white noise, upper 
cut-off frequency of /„ = 1024 Hz and observation time of one 
sidereal day. From the plot it is evident that the singular val- 
ues are negligible after ~ 1000. These results agree with the 
numerical and SPA results which give the size of the central 
spot ~ 0.1 radian. The dotted line shows the singular values 
for the same detector pair, but with LIGO-I noise PSD and 
upper cut-off frequency of f u — 512 Hz. Clearly, the latter 
curve indicates a loss of angular resolution due to relatively 
poorer higher frequency response of the detectors. 



The width of the beam estimated above seen to be con- 
sistent with the 'number of degrees of freedom' present in 
the kernel (beam). A widely used method that identifies 
the linearly independent modes in a linear transformation 
is the singular value decomposition (SVD) [25]. The de- 
composition identifies linear combinations of modes that 
have almost zero eigenvalues - the null subspace. It then 
projects out the solution orthogonal to the null subspace 
which spans the true degrees of freedom. The singular 
values of the kernel for the LIGO detectors at Hanford 
and Livingston for white noise with upper cut-off fre- 
quency /„ = 1024 Hz are plotted in Figure 5 (solid line). 
The figure shows that the eigenvalues become essentially 
negligible after ~ 1000 implying that this is the number 
of degrees of freedom in the kernel. The numerical and 
the SPA analysis shows that the size of the central spot 
or the resolution is ~ 0.1 radian, which means that there 
are 47r/(0.1) 2 ~ 1000 independent patches in the sky. So 
the SVD results are consistent with the size of the beam 
obtained by numerical and theoretical methods. 

In practice, however, the higher frequency response of 



IV. THE MAXIMUM LIKELIHOOD 
DECONVOLUTION 

A. Unpolarised Background and Single Baseline 

We first consider the simpler case of detecting and de- 
convolving the signal from an unpolarised SGWB using 
one pair of detectors. We later indicate in the subsections 
that follow how this method can be extended to the more 
general cases of SGWBs and detector baselines. 

The observed data construct, S(Cl), consists of a signal 
and additive noise, namely, 



S((l) = a(«) + n(«) 



(4.1) 



The first term on the rhs, the expectation of the dirty 
map given in Eq. (3.37), is a convolution of the true power 
in SGWB Va(Q) in the two polarizations arriving from a 
direction f2 in the sky with corresponding beam response 
functions B A (fl,fl'). Note that, though the definitions 
of the above quantities involve complex Fourier trans- 
forms, these quantities are all real owing to the fact that 
they were originally derived from real time series, so that, 
SJ(t; -/) = 3>(i; /), 7*(«, t, -/) = 7 (0, t, f) and so on. 

In this section, we construct the Maximum Likelihood 
(ML) estimator for the angular power distribution Va{&) 
given the measured data S(fl). For simplicity and clarity 
of presentation, we first limit the analysis to the simple 
case where both the polarizations follow the same angular 
power distribution V(tl). This simplifies the form of the 
construct to 



s(n) = / dn' [B+(n,n')+B x (n,n')]v(n')- 

Js 2 



• n(n). 

(4.2) 

In practice, the sky is divided into finite number of 
pixels. Then, the observed data vector is denoted by S, 
whose component Si := S(fli) is the signal measured at 
the i th pixel. We similarly define the vectors "P and n, 
with components 'P(fij) and n(fli), respectively. In this 
notation the convolution leads to a set of linear algebraic 
equations, 



B V 



(4.3) 



where B is the known [51] beam matrix, expressed as 
Bij := B + (tti, ttj) + B x (tti, ttj). In the weak-signal ap- 
proximation, the variances of the signal-noise cross terms 
^i 2(^1 /) ^2,1 it] f) are much smaller than the variance of 
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the noise-noise cross term n\ (t; /) n,2 (t; f) . So the ob- 
served pixel noise is strongly dominated by the noise- 
noise term and can be written as, 



n(«i) 
dt 



J diA^O,^ 
dfnl(t;f)n 2 (t;f) 



x (4.4) 
Q(fli,t,f;H) t 

A(t,no ' 



the cross-terms h* 2 {t; f) ri2,i(t; /) have been dropped. 
The pixel noise is a sum of a large number of zero mean 
random numbers, where none of the addend strongly 
dominate (statistically) over the others. Hence, follow- 
ing the generalized central limit theorem [24], one can 
argue that the pixel noise tends to be zero mean Gaus- 
sian. If this argument is used to calculate the variance 
of n(Cl), we consistently get the same result as expected 
from Eq. (2.17). After a routine but fairly involved alge- 
bra, the pixel-to-pixel noise covariance matrix, N = Nij 
of the dirty map turns out to be: 



(m rij) 



At 



dt A _1 (i, flj) 



Bij. (4.5) 



The deconvolution problem is, of course, a very stan- 
dard problem in many areas of science. In particular, 
Eq. (4.3) is identical in structure to the set of equations 
that arise in the map-making stage of CMB experiments. 
The temperature anisotropy ATi in a direction is in- 
ferred from time-stream data, d t using a linear model 
d t = '^2 i A t iATi + n t . The convolution kernel A t i that 
relates the time domain to the pixel domain is determined 
by the pointing or scan strategy as well as the beam re- 
sponse function of the antenna. The noise n t is (assumed 
to be) Gaussian and described by the noise covariance 
matrix N tt i. As described below, a ML solution for the 
sky map ATi is readily obtained in this linear model un- 
der the assumption that the noise is Gaussian. We adapt 
this technique to solve our problem since it has been ap- 
plied with great success in the CMB field and there exists 
an extensive literature [26] and also public domain pack- 
age for implementing it numerically [27, 28]. However, 
it should be noted that the problem differs in two im- 
portant aspects. First, in our case there is the simplicity 
that the kernel connects two vectors which are both in 
sky pixel space. This implies the kernel is a square ma- 
trix for the case of single baseline and unpolarized back- 
ground. Second, in our the case the statistics of the noise 
is potentially non-trivial; the noise in this situation is a 
complex object built by integrating the product of two 
random variables corresponding to the noise streams in 
each detector. The Gaussianity of the noise has to arise 
from the generalized central limit theorem [24]. 

We proceed assuming that the joint probability distri- 
bution of the elements of n is a multivariate Gaussian dis- 
tribution [29] given by the probability distribution func- 
tion: 



<P(n) 



1 



(27r)*W 2 ||N|| V2 
1 



exp 



-n 1 - NT 



( 27r )iW2 



x exp 



(n^rST 1 -n+Tr[lnN]) 



(4.6) 



where iV P i X is the total number of pixels and N := (nn T ) 
is the known noise covariance matrix [see Eq. (4.5)]. 
Thus, given a signal vector TP, the probability of observ- 
ing radiometer output S is [27], 



<P(S|7>) - (2tt) 



exp 



-Afp.x/2 



; ((S - B • V) T ■ INT 1 ■ (S - B • V) 



Tr[lnN]) 



(4.7) 



Solving for d^(S\V)/dV = 0, using the fact that N is 
symmetric and positive definite, it is straightforward to 
show that the above probability is maximum when 



V = (B T N _1 B) * B^'IST 1 • S 



T Tvr-1 



(4.8) 



which is, therefore, the well-known result for the desired 
ML estimator of the true sky map of SGWB anisotropy. 
The deconvolved map will also have pixel noise, given by 

n := V V 

= (B T N _1 B) -1 B T N" 1 • (B V + n)- V 



B r N _1 B) " B' J 'N 



T Tvr-1 



n 



(4.9) 



and the pixel to pixel noise covariance matrix of the ML 
map is obtained as 



S := 



(nn 



= (B 



T N- 1 B) 



(4.10) 



Therefore, to obtain the ML map estimate one has to first 
compute the inverse of the pixel to pixel noise covariance 
matrix XT 1 = B T N X B. The ML map is then obtained 
as a solution to the linear algebraic equations, 



1T X V = B N _1 • s. 



(4.11) 



In solving this equation, we have a choice of either 
using one of a number of direct methods or one of the 
iterative methods. In the low resolution regime, the ma- 
trix inversion looks feasible with reasonable accuracy. 
However, in general, direct methods are computation- 
ally more expensive and iterative methods are preferred, 
provided their convergence to the solution is rapid. For 
sparse matrices the iterative Conjugate Gradient (CG) 
method is well suited for this inversion [25, 30] problem. 
The CG method solves linear systems with symmetric 
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positive definite matrix and have been found to be more 
advantageous compared to other iterative methods such 
as the Jacobi method [26, 31]. Starting with a guess solu- 
tion, the convergence of the method can be often greatly 
improved by 'preconditioning' the system of equations, 
i.e., multiplying both sides with a suitable matrix (say, 
the inverse of diagonal elements of the matrix.). Our 
choice is also motivated by the fact that the CG method 
has been successfully implemented for map making in 
CMB experiments [26, 27, 32, 33, 34]. 

The above method can be extended to include multiple 
baselines and also to estimate power in each polarisation 
component. We discuss these extensions in the following 
subsections. 



with X representing the matrices S, B, and n. Note 
that S, n are now 1 x N p - lx Nh vectors and B is a iVpj x x 
Npi x Nb matrix, while T 3 (the true SGWB sky) remains 
unchanged. This is similar to CMB experiments where 
each pixel is visited by the detector several times. In the 
multi-baseline GW radiometer case each pixel is visited 
by different baselines, and, unlike a CMB experiment, 
each pixel is visited an equal number of times. In this 
case too the Maximum Likelihood estimate [52] and the 
pixel to pixel noise covariance matrix of the ML map are 
given by: 



B. Multiple Baselines 



The above analysis can be extended to a set of 
GW radiometer baselines. Let be the observed map 
by the i th baseline with beam matrix B^ and observed 
noise n^. Then Eq. (4.3) can still be written as: 



SB' 1ST 



B T N _1 B, (4.13) 



BP 



where 



X 



X (2) 1 



V xw / 



(4.12) 



but the noise covariance matrix N := (nn T ) of the raw 
sky map has to be modified. Let n^tt) be the pixel noise 
from the radiometer baseline i with detectors I and I' 
and we follow the same convention for Ai and Qi. Then, 



(n i (n 1 )n J (^ 2 )) = 



dt A^ftOi) 



dh / dt 2 



d/i 



d/ 2 



dtxj\t,n 2 ) 



d/i / d/^ A *(/i-/i)<M/2-/2) x 



{n I {ti,fi)nj{t 2 ,f2)nii {h, / 1 )n J /(f 2 , / 2 ))- 



A 3 -(ta,n 2 ) 



(4.14) 



If i and j denote the same baseline we get back the previ- 
ous result [Eq. (4.5)]. However, if i and j denote different 
baselines, at least one of the two detector pairs will be 
different (i.e. cither / ^ J or 7' ^ J'), so in that case 
(rii(tti)nj((l 2 )) = 0. Hence, one can write 



(4.15) 



Thus, the matrix N for a network of baselines will be a 
block diagonal matrix, where each diagonal block is 
the noise covariance matrix for the corresponding base- 



line, N, 



(<) 

kk' 



(nj(f2fe)rij(f2fe')) , provided in Eq. (4.5). 



The above algebra suggests that it is fairly straight for- 
ward to combine the observations made by multiple base- 
lines for the estimation of the true SGWB angular power 
distribution. In [17, 35] it was shown that the optimal 
way to search for a isotropic SGWB using a network of 
detectors (with uncorrelated noise) is to linearly combine 
correlations from pairs of detectors, instead of comput- 
ing higher order correlations using data from more than 
two detectors. We can extrapolate the same logic to the 
directed search and argue that the procedure described 
above to combine data from a network of detectors is also 
optimal. 
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The search for a SGWB using a network of detectors is 
becoming progressively relevant as other kilometer scale 
detectors, namely Virgo, GEO and LCGT, are expected 
to reach their initial goal sensitivity in the next few years. 
A network of detectors can enhance the directed search 
in many ways. The resolution of a radiometer is propor- 
tional to the length of the baseline. Inclusion of a detec- 
tor at a distance like Virgo, which is further away of the 
LIGO detectors than the mutual separation of the LIGO 
sites, will clearly increase the highest baseline separation 
and hence the resolution of the baseline. However, more 
important enhancement would be realized due to better 
coverage of the sky. An analogy with radio astronomy us- 
ing an array of antennas may be appropriate to mention 
in this context [53]. As the earth rotates, the projections 
of the radio antenna baselines on the plane perpendicu- 
lar to the source direction sample different points on the 
two dimensional Fourier plane (commonly known as the 
u-v plane). The sampled Fourier plane is then inverse 
transformed to generate the image. While the highest 
resolution of the network is limited by the projection of 
maximum antenna separation, addition of more anten- 
nas to the network ensures better sampling of the u-v 
plane reducing the side lobes, thereby producing a more 
faithful image of the sky. Detailed introduction to the 
basic principles of earth rotation aperture synthesis can 
be found in most of the standard texts on radio astron- 
omy, e.g., [36]. In GW radiometry with a network of 
detectors we expect that a similar scenario will arise - 
better coverage of the sky should be possible due to dif- 
ferent orientations of the baselines with respect to the 
source. Moreover, since the true power distribution will 
be estimated from an over-constrained set of equations, 
the error in the estimated quantities will be reduced. In 
addition, a radiometer search can benefit from certain 
technical advantages that a network of detectors can pro- 
vide: A detector at a third site joining the LIGO detec- 
tors will boost the "single-baseline integration time" , i.e., 
the single-baseline duty-cycle in a three-site network will 
be at least as good as, but will likely be better than, 
that in a two-site network. Also, owing to common in- 
strumental noise sources in the LIGO detectors, certain 
frequency bands are currently notch-filtered in comput- 
ing the cross-correlation statistic. Some of these noise 
sources are known not to affect Virgo and, therefore, the 
LIGO- Virgo cross-correlation statistics. For example, the 
power-line noise affects the LIGO detectors at the har- 
monics of 60Hz, whereas it affects the Virgo detector at 
the harmonics of 50Hz. Therefore, a radiometer search 
that benefits from the LIGO-Virgo baseline's contribu- 
tion will probe the presence of astrophysical signals over 
a larger set of frequencies than one limited to the baseline 
consisting of the LIGO pair of detectors. 



C. Polarization Map 

We may also choose to extract power from different 
polarizations separately. The discrete convolution equa- 
tion. 



B 



B, 



(4.16) 



can also be expressed by Eq. (4.3) 
S = B V + 

where 



B := B, B> 



V := 



(4.17) 



Again in this case the ML estimator of sky map can be 
expressed by Eq. (4.13): 



V = SB' N 



T-vr-1 



XT 1 := B T N _1 B 



This case can also be generalized for a network of de- 
tectors by retaining the same definition of V [Eq. (4.17)], 
but redefining the beam matrix as: 



B 



(i) 
+ 

(2) 



(1) 

X 

(2) 



B (N b ) B (JV b ) 



(4.18) 



and, again, using the same ML estimation formula given 
in Eq. (4.13). 



V. IMPLEMENTATION AND NUMERICAL 
RESULTS 

In this work we have numerically implemented the 
maximum likelihood estimation algorithm on simulated 
data using the MATLAB® software package [37] to esti- 
mate the "true" (unpolarized) SGWB sky observed with 
a single baseline ground based GW radiometer. The de- 
tails of the computation scheme, the numerical deconvo- 
lution algorithm, the simulated data and the deconvolved 
maps are presented in this section. 



A. Preparation of Simulated Dirty Maps 

The data are simulated in the frequency domain for 
each time segment. Since the noise of ground based in- 
terferometric detectors is very high at frequencies greater 
than few 100 Hz and the computation cost increases with 
the number of frequency bins, we use an upper cut-off 
frequency of f u — 512 Hz and bin width of A/ = 2 Hz 
for testing of the algorithm. The length of each time 
segment is chosen as At = 192 sec and the total inte- 
gration time is T = 86400 sec [54]. The sky is pixelized 
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using the Hierarchical Equal Area isoLatitude Pixeliza- 
tion (HEALPix) [55] scheme, which divides the 2-sphere 
(S ) in 12 n;? ide pixels, where n s ^ c is an integer power 
of 2. Since the radiometer beam-width is greater than 
~ 6°, we chose n s id c — 16, which corresponds to a pixel 
width of - 3° and a total of 3072 pixels. The HEALPix 
scheme also allows fast spherical harmonic transform on a 
sphere, which may become useful for more advanced anal- 
ysis in future. Note that, the algorithm is independent 
of the pixelization scheme, other equal area pixelization 
schemes can also be used in the analysis. 

We generate the detector noise nj(t; /) using a Gaus- 
sian pseudo random number generator for each time seg- 
ment. The noise is colored using the (one sided) noise 
PSD Pi(t; /) of the corresponding detector according to 
Eq. (3.15). MATLAB® software's pseudo random num- 
ber generator randn can generate very long sequences of 
random numbers, so we relied on that routine for simulat- 
ing detector noise. For each of T/At = 86400/192 = 450 
time segments, we generated a complex random sequence 
(that is, two real random sequences) of f u /Af = 512/2 = 
256 real numbers. The total number of random numbers, 
2(T/At)(f u /Af) = 225, 000, is much less than the period 
of randn, which is 2 1492 > 10 449 [38]. 

Signal is also generated directly in the frequency do- 
main. However, the GW strain in each detector, hi(t; /), 
is not generated independently; rather the product of 
the strains in the detectors, h*(t; f) hi{t\ /), is gener- 
ated directly using the statistical properties of the strain 
correlation described in subsection III A, in particular, 
Eq. (3.10). We may write the product of the strains as a 
sum of its expectation value and statistical fluctuations: 

K(tJ)h 2 (tJ) = (hl(tj)h 2 (tj)) + fluctuations. 

(5.1) 

Since our main aim is to generate, 

T, (t , /) s 2 (t , /) = [hi (t , /) + n\ (t , /)] [h 2 (t , /) + n 2 (t , /)] 

(5.2) 

and since statistically the variation in the signal terms are 
much weaker than the zero mean uncorrelated detector 
noise terms, we may simply drop the signal "fluctuations" 
term from Eq. (5.1) - that is, we may approximate the 
product of the detector outputs using the formula[56]: 

s{{tj)s 2 {t,f) = (hl(t,f)h 2 (tj)) + nl(t,f)n 2 (tj). 

(5.3) 

For all the cases considered in this paper we have used 
flat source PSDs, i.e., H(f) — constant. 



In this analysis we assume the sky to be a collection of 
uncorrelated point sources of different strengths placed at 
every pixel. Moreover, the numerical analysis has been 
restricted to the case of equal power in each polarization. 
So the (injected) true sky is constructed by putting 

Ptrue(n) = J2 V kS(n-n k ), (5.4) 

k 

where Vk is the strength of the point source placed at 
pixel k, located in the direction flj. (in order to inject 
only one point source at pixel ko, we make all the Vk = 
except for k = ko). In this set up, the expression for the 
overlap reduction function [Eq. (3.7)] for the true SGWB 
strain becomes 

j(t,f;An,dV A ) = ^rfft^l^e 2 """'^"^, 

k 

(5.5) 

where we use our usual notation for T(Ct,t) defined by 
Eq. (3.34). 

We substitute the above in Eq. (3.10) and inject that 
simulated signal in noise using Eq. (5.3) to generate prod- 
ucts of outputs from two detectors. In order to pre- 
serve the reality of time series data, the products of sig- 
nals are generated only for positive frequencies and set- 
ting the negative frequencies equal to the complex con- 
jugates of their positive frequency counterparts, that is, 

»/(*,-/) =3}(*,/). 

The radiometer analysis is then run on the simulated 

data by placing filters Q(ttk,t, /; H) at each pixel k to 

generate the dirty maps. 

B. Deconvolution: Clean Maps 

Any deconvolution routine requires the beam function 
at all the points of interest on the sky. For a GW ra- 
diometer, the filters (and hence the beam functions) are 
dependent on the data set itself. So, if the beam function 
is calculated for each sky pixel, apparently the computa- 
tional cost should go up by a factor of number of pixels 
(^Vpix ~ 3000) times the cost to make one sky map of 
beam for a given pointing direction. However, we can 
use algebraic tricks to make this method computation- 
ally viable. A possible way to implement this method for 
the simple case of one baseline and equal power in each 
polarization is demonstrated below. 

The beam and noise covariance matrices are given by: 
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(5.6) 
(5.7) 



r 



where 
G(tJ) 



F+(ti, t)F+(n, t) + f* (n, t)F* (n, t) , 



ff 2 (/)/[Pi(«,/)P 2 (t,/)] 



(5.8) 



Now one can see that the beam matrix is a summation 
of parts which depend on either fi^ or f2j. So, it is pos- 
sible to precompute the arrays \{tl,t), T{£l,t), G(t,f) 
and Ax(t), only once, and then use them efficiently to 
evaluate the whole B matrix. Since each of these arrays 
are functions of any two of the three variables tt,t 7 f, the 
memory size required for each two dimensional array will 
not be too large. Once the arrays are computed, the next 
step is to obtain each clement of the beam matrix sepa- 
rately. Each element requires the evaluation of two sum 
loops (like matrix multiplication) involving one exponen- 
tial. The symmetry of the beam matrix (without the 
normalization constant) can be utilized here to reduce 
the computation cost by a factor close to 2. To summa- 
rize, a significant amount of CPU time required to make 
a skymap of the beam for one pointing direction could 
be utilized to make the beam maps for all other pointing 
directions. This is true for the noise covariance matrix as 
well, which, in this case, is proportional to the (unnor- 
malized) beam matrix. A Fast Fourier Transform (FFT) 
with interpolation trick [42] and assumption of station- 
arity of noise for deconvolution ( not for making the dirty 
map, where non-stationarity will be accounted for) can 
also be incorporated in future to reduce the CPU time. 

In this simple case, since the beam matrix is a square 
matrix, so that, £ = (B T N- 1 B)" 1 = B- 1 N(B- 1 ) T , 
the estimated map given by Eq. (4.13) is just the least 
square solution: 



B 



(5.9) 



Even for the general cases of multiple baselines and polar- 
ized background, the estimation equation takes the above 
simple form [see section IV]. 

The first task was to compute the beam matrix, which 
happens to be the computationally most intensive task. 
A typical beam matrix for the LIGO baseline using 192 
HEALPix pixels is shown in Figure 6. By construction 
Bkk = 1 and |-Bfcfc'| < 1 for k ^ k' ', hence the matrix 
is diagonal dominated. The "stripes" in the matrix are 
related to the pixclization scheme. The beam is stronger 
if the pixels are closer to the pointing direction and it 
weakens as the distance between the pixels and pointing 



direction increases. In other words, the pixels closer to a 
point source will have stronger contamination from the 
point source. However, since we have used a isoLatitude 
pixelization scheme, the indices of two neighboring pixels 
at different latitudes differ by the total number of pixels 
on that latitude. This fact is reflected in the plot of the 
beam matrix - the matrix is sparse with certain "peri- 
odic" behavior which produces the stripes in the plot. 
The matrix becomes even more sparse for finer resolu- 
tions as greater number of isoLatitude pixel rings pass 
through the core of the beam. Making a legible plot of 
the beam matrix for higher resolution is difficult, so the 
plot presented here is restricted to lower resolution - 192 
pixels instead of 3072. 
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FIG. 6: A typical Beam matrix for the LIGO Hanford- 
Livingston baseline at a low resolution (192 pixels) is shown 
in this figure. Each row of the matrix is the beam response 
function for the pointing direction that corresponds to the row 
index. The matrix is diagonal dominated as the radiometer 
receives maximum contribution from the pointing direction. 
The stripes are related to the isoLatitude pixelization scheme 
- the indices of the neighboring pixels at different latitudes 
differ by the total number of pixels on that latitude. The 
possibility of making the beam matrix more diagonal domi- 
nated using a nested pixelization scheme, where the indices 
of the neighboring pixels are close, is being explored. 

Since the sparseness of the beam matrix depends on the 
pixelization scheme, it may be possible to make the beam 
matrix significantly diagonal by using a nested pixeliza- 
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tion scheme, where the indices of the neighboring pixels 
are close. This possibility is being explored. 

Sparse matrices are computationally easier to invert; 
however, the stability of inversion of such matrices is 
a numerical challenge. Therefore, as mentioned in sec- 
tion IV, instead of evaluating P — B 1 • S [Eq. (5.9)], we 
choose to algebraically solve for P from the set of linear 
equations 



BP 



(5.10) 



with the same number of unknowns as the number of 
equations (i.e., the system is not under or over con- 
strained). The above equation can be cast into a linear 
system with symmetric kernel, 
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(5.11) 



by introducing two new quantities, 
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(5.12) 

Comparing with Eq. (5.6) one can clearly see that b is 
a symmetric matrix. A symmetric kernel is always pre- 
ferred by most of the algorithms for solving numerical 
linear equations. Moreover, it is possible to show from 
the algebra presented in subsection IV B that, to incorpo- 
rate observations from multiple baselines, one can simply 
replace s and b by the sum of those respective quantities 
over all the baselines and solve Eq. (5.11) for P (which is, 
of course, independent of any baseline). Thus the exten- 
sion of this analysis to incorporate a network of detectors 
becomes straight forward with these new quantities. 

Following the discussion presented in section IV, we 
use a Conjugate Gradient (CG) iterative technique to 
solve the above set of linear equations. CG algorithm 
solves a set of linear equations A • x = b, where A 
is a square matrix and x, b are vectors, by minimizing 
the quadratic form | x • A • x — b • x. We use the 
minimum residual method, which efficiently utilizes the 
fact that b is symmetric and does not require b to be 
positive definite. The minimum residual method aims to 
minimize the residual |A • x — b| 2 itself, instead of the 
quadratic form |x ■ A ■ x — b ■ x. Further details on CG 
methods can be found in standard literature, e. g., [39]. 

The clean maps also contain pixel noise - partly due 
to the random noise present in the data and partly due 
to the numerical errors introduced at each stage of the 
pipeline, mainly during the process of deconvolution. 
There are pixels in the deconvolved map, which have neg- 
ative values, even though the injected map is positive. To 
reduce the noise in the clean maps, we introduce an ad- 
ditional step: We compute the root-mean-square (RMS) 
noise "cr" in a map when there is no injected source. Then 
in the clean map (with source) we mask, that is, set to 
zero, all the pixels that have values less than a thresh- 
old of few a. The number of iterations for deconvolution 
and the threshold for masking can be adjusted according 



to the tolerable levels of false alarm and false dismissal 
probabilities. 

The above can be easily extended to handle real data 
where we have no control on the injections. One can cal- 
culate the equivalent of RMS noise for no injection by 
shifting the data streams from different detectors by a 
large time lag (much smaller than the segment duration) , 
say, 1 sec, that corresponds to distances much greater 
than the earthly distances, so that, true GW signals are 
not added coherently. To be more careful, one can per- 
form this exercise for a few large time shifts and confirm 
that the noise levels are not significantly different for dif- 
ferent shifts. 

It is, however, not so straight forward to measure the 
quality of deconvolution. The signal-to-noise ratio (SNR) 
of individual pixels do not carry enough significance, as 
the neighboring pixels are highly correlated. It is also 
difficult to define a quantity that can take into account 
the pixel-to-pixel noise covariance due to the difficulty in 
inverting the beam matrix. In this paper, we use a rather 
simplistic measure to quantify the quality of deconvolu- 
tion, which is often used in image processing to measure 
the reconstruction error. We use a quantity known as 
the "Normalized Mean Square Error" (NMSE) [40], ex- 
pressed in terms of the injected P and the estimated P 
maps as 



NMSE := 



IP -PI 



(5.13) 



Obviously, lower the NMSE, better the reconstruction. 

The whole analysis was tested for different kinds of 
injected maps consisting of localized sources and diffuse 
sources. In all these cases each pixel k of the injected 
map was assigned a value Vk between and 1 with a 
source PSD H(f) = 5 x 10~ 47 /Hz. This means that, if 
a pixel of a test map has strength 1, the standard devi- 
ation of the Fourier transform of stochastic GW coming 
from that pixel is y/H(f) - 7 x 10~ 24 /\/i?z. This stan- 
dard deviation is about one third of the standard devi- 
ation of Fourier transform of noise at the most sensitive 
frequency band of the LIGO-I detectors which is about 
2 — 3 x 10~ 23 j\jHz. To our knowledge, the strength 
of anisotropic astrophysical GW background has not so 
far been predicted theoretically. However, if we try to 
extend the results from the all-sky averaged (isotropic) 
astrophysical background [41] to have a crude estimate 
of the strength of the anisotropic background, it turns 
out that, the PSD of the anisotropic astrophysical back- 
ground in the universe is weaker than the H(f) we have 
injected by roughly a few orders of magnitude. In the 
present work we have used an observation time of one 
day to demonstrate the method. With longer observa- 
tion times and employing several baselines comprising of 
the upcoming advanced (more sensitive) detectors, the 
difference between the expected background and the de- 
tectable background would diminish or altogether disap- 
pear. 



20 




0.2 0.4 0.6 0.8 1 

(a) Injected map 
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(b) Dirty map - made from simulated data using the 
radiometer analysis 
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(c) Clean map - obtained by deconvolution of the dirty 
map using 15 CG iterations (NMSE = 1.22) 
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(d) Masked clean map - all the pixels below a threshold 
of 5<r were set to zero (NMSE = 0.64) 



. 7: Illustration of deconvolution for localized sources: A 
pixel wide localized source was injected near the location of 
Virgo cluster - a potential source of SGWB. 



We first inject a 4-pixel wide localized source near the 
Virgo cluster, a potential point source of SGWB, as il- 
lustrated in Figure 7. Figure 7(a) shows the injected 
map. Figure 7(b) shows the dirty map made from simu- 
lated data. Figure 7(c) shows the clean map, obtained by 
deconvolving the above dirty map with the beam using 
15 conjugate gradient iterations and Figure 7(d) shows 
the same clean map masked using a 5<7 threshold. It is 
evident that deconvolution has successfully localized the 
source in a relatively smaller area as compared to the 
dirty map. Still, one should note that, the deconvolution 
routines do not perform well when the injected source is 
like a delta function. This causes high NMSE, 1.22 for 
the unmasked and 0.64 for the masked clean maps, and 
significant loss of the peak strength of the reconstructed 
point source, as indicated in Figure 7. Moreover, increas- 
ing the number of iterations beyond a certain level actu- 
ally deteriorates the quality of deconvolution due to noise 
amplification, and this level is dependent on the kind of 
source one is searching for. In this basic analysis we have 
used 15 iterations to search for localized sources and 40 
iterations to search for broad sources, which offer reason- 
ably clean deconvolution and comparatively low NMSE. 
Introduction of a minimum error criterion [43] to termi- 
nate the iteration process is being considered. Several 
other deconvolution algorithms are being explored in or- 
der to identify the one which is best suited for the GW 
radiometer analysis. 

Next, we inject two kinds of diffuse sources, viz., one 
that is nearly equatorial and another that is distributed 
across multiple declinations, as illustrated in Figure 8. 
We injected modified (using FTOOLS[57]) galactic fore- 
ground seen in CMB temperature anisotropy measure- 
ments as our test patterns for the diffuse SGWB sources. 
The left panels of Figure 8 correspond to a modified 
form of the temperature anisotropy map measured by the 
WMAP satellite [44]. We emphasise that the sky looks 
different in barycentric coordinates, similar to what is 
shown in the top- right panel of Figure 8. We omit the 
coordinate transformation step intentionally in order to 
get a diffuse equatorial source. The right panels of Fig- 
ure 8 correspond to a modified form of the temperature 
anisotropy map in the barycentric coordinates generated 
by the Planck Simulator [46] . One of the main modifi- 
cations applied to both of these maps was to mask the 
brightest part of the galaxy. This step reveals more struc- 
tures in the maps, which is useful for testing a deconvolu- 
tion algorithm. Figure 8(a) shows the injected toy maps. 
Figure 8(b) shows the dirty maps obtained by the ra- 
diometer analysis. One can see that the dirty maps have 
lost all the fine structures present in the injected maps. 
Furthermore, they show certain features that were not 
even present in the injected map. Also, the pixel val- 
ues in the dirty maps are spread over a range consist- 
ing of large positive and almost equally negative values. 
Figure 8(c) shows the clean maps recovered by 40 CG 
iterations. Clearly, many of the features of the injected 
maps have been recovered in the clean maps, which is also 
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(c) Clean maps - obtained by deconvolution of the dirty maps using 40 CG iterations. Clearly, the structures and positivity 
of the injected maps, which were lost in the dirty maps, have been restored quite significantly. 

(left panel: NMSE = 0.33; right panel: NMSE = 0.22) 
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(d) Masked clean maps - obtained by putting a 4cr threshold on the clean maps 
(left panel: NMSE = 0.36; right panel: NMSE = 0.33) 



FIG. 8: Illustration of deconvolution for broad sources: Maps similar to CMB temperature anisotropy sky with the galactic 
foreground were injected as toy maps. The left panels correspond to a map measured by the WMAP satellite [44], as seen in 
the Galactic coordinates, as a toy model for an equatorial source and the right panels correspond to a map generated by the 
Planck Simulator [46], as a toy model for a multi-declination extended source. 
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evident from the lower values of NMSE, 0.33 and 0.22 
respectively. Also, the positivity of the estimated map 
has been vastly improved - the pixel values of the clean 
maps lie mostly on the positive side, as one should ex- 
pect. Finally, Figure 8(d) shows the masked clean maps 
obtained by using a 4cr threshold. Though the masked 
maps give better visual impression, masking can actually 
discard several pixels which have weak sources, thereby 
increasing the NMSE. In Figure 8(d), for example, mask- 
ing increases NMSE to 0.36 and 0.33 respectively, though 
the masked maps look more similar to the injected maps 
shown in Figure 8(a), than the unmasked clean maps in 
Figure 8(c). 

VI. CONCLUSION 

The stochastic astrophysical GW background is likely 
to be dominated by sources in the nearby anisotropic uni- 
verse, so the detection of localized sources is more favor- 
able than the all-sky-averaged search. Making a skymap 
of the SGWB sky has been a long standing ambition 
of stochastic GW research. Different analysis methods 
have been proposed to create skymaps by measuring the 
first few spherical harmonic multipoles of the sky. Here 
we have presented a direct approach of directed GW ra- 
diometer analysis. In this approach, the whole sky is 
decomposed in a discrete set of pixels and the contri- 
bution from each pixel is measured separately by cor- 
relating phase shifted detector outputs to generate the 
whole skymap, which is a clear application of earth ro- 
tation aperture synthesis. Specifically, for the AGWB 
detection statistic, we have defined a correlation statistic 
with a directed optimal filter that targets a fixed point 
in the sky by adjusting the time-delay across a baseline 
to track its rotation with the Earth. This statistic how- 
ever provides us with a dirty map of the sky which we 
numerically deconvolve to obtain the true skymap. For 
this purpose, we have employed the conjugate gradient 
method. We numerically implement the deconvolution 
on simulated unpolarised GW sky maps obtained with 
the LIGO detector baseline. The success of this method 
is demonstrated by the recovery of simulated source dis- 
tributions, namely, (i) of a point source, (ii) of a diffuse 
source in the equatorial plane, and (iii) of a diffuse source 
off the equatorial plane. 

This work needs to be implemented on other base- 
lines of the upcoming/future network of detectors such 
as LIGO- VIRGO, LIGO-LCGT, LISA etc. The outline 
of the analysis has been presented here. However, fur- 
ther detailed analysis and the implementation is a future 
goal. Even for the single baseline, SPA analysis shows 
that, perhaps a more efficient method of deconvolution, 
yielding better accuracy and convergence, lower com- 
putational costs and convenience of application may be 
possible using more sophisticated analysis, for example, 
one involving basis functions. The Maximum Likelihood 



framework presented here is, in fact, independent of any 
particular choice of basis. So once a suitable basis is cho- 
sen, rest of the analysis can be applied without requiring 
any major change. It may also be possible to deconvolve 
only a patch of the sky using a similar method [45] . 

The work presented in this paper should also benefit 
two other searches. First, since the long-duration integra- 
tion of the data will essentially comprise a sum over short 
stretches, a large signal in a short stretch will constitute a 
candidate for a transient or burst (short-duration) event. 
Unlike the coincidence search being currently conducted 
for such events, our work will combine coherently the 
outputs of several detectors and, thus, improve their de- 
tectability. Second, the long-duration integration of the 
data should be able to find gravitational wave signals 
from modelled sources, such as pulsars. Although our 
method is optimal for searching unmodelled sources, it is 
not so for pulsars, the signals from which can be matched 
filtered. The problem with the latter method is that ow- 
ing to the very large parameter space volume, an all-sky, 
all-frequency search for pulsars with matched filtering is 
not computationally viable. Our proposed method is not 
handicapped by this problem since it does not use the 
intrinsic source parameters for the search; rather it uses 
the data from one detector to 'filter' that from others 
in the network, after appropriately time-shifting them 
and weighting them with the respective antenna patterns. 
Thus, our method can be used as the first step in a two- 
step hierarchical search for pulsars, where triggers from 
our method are followed up with matched filtering. 
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